




#ifndef _C_MATRIX_H_
#define _C_MATRIX_H_




#include "sl_complex.h"



/// <summary>
/// C语言接口的向量结构体
/// </summary>
typedef struct _vector__
{
	/// <summary>
	/// 向量的数据
	/// </summary>
	real_t* data;
	/// <summary>
	/// 向量元素个数
	/// </summary>
	int size;
}vec;


/// <summary>
/// C语言接口的向量结构体
/// </summary>
typedef struct _vector_int_
{
	/// <summary>
	/// 向量的数据
	/// </summary>
	int* data;
	/// <summary>
	/// 向量元素个数
	/// </summary>
	int size;
}ivec;



/// <summary>
/// 矩阵的C语言结构体接口，矩阵元素全部是real_t型变量。
/// 注意，只有矩阵乘法有函数内并行加速的可行性。其他求逆，QR分解，奇异值分解等，只能在函数外部进行并行加速。
/// </summary>
typedef struct _matrix__
{
	/// <summary>
	/// 矩阵的二维数组指针。是实数类型的变量。
	/// </summary>
	real_t** data;
	/// <summary>
	/// 矩阵的行数和列数
	/// </summary>
	int rows, cols;
}mat;




//矢量内存分配。 v：矢量对象（这里必须用引用，因为要修改指针指向）。size：矢量大小。
void vec_allocate(vec* v, int size);
//回收矢量内存
void vec_deallocate(vec* v);


//矢量内存分配。 v：矢量对象（这里必须用引用，因为要修改指针指向）。size：矢量大小。
void ivec_allocate(ivec* v, int size);
//回收矢量内存
void ivec_deallocate(ivec* v);


//内存分配。 
void mat_allocate(mat* m, int rows, int cols);
//回收内存
void mat_deallocate(mat* m);
//

/// <summary>
/// 向量打印函数
/// </summary>
/// <param name="v">打印的向量</param>
void vec_print(vec v);


/// <summary>
/// 向量打印函数
/// </summary>
/// <param name="v">打印的向量</param>
void ivec_print(ivec v);

/// <summary>
/// 控制台打印矩阵
/// </summary>
/// <param name="m">矩阵</param>
void mat_print(mat m);


/// <summary>
/// 向量随机化
/// </summary>
/// <param name="v">待随机化的向量</param>
/// <param name="seed">随机数种子</param>
/// <param name="min">随机数的最小值</param>
/// <param name="max">随机数的最大值</param>
void vec_randomize(vec v, int seed, real_t min, real_t max);


/// <summary>
/// 向量模长归一化
/// </summary>
/// <param name="v">归一化的向量</param>
void vec_normalize(vec v);

/// <summary>
/// 向量初始化
/// </summary>
/// <param name="v">初始化的向量</param>
/// <param name="r">初始化的值</param>
void vec_init(vec v, real_t r);

/// <summary>
/// 实现向量的复制。
/// </summary>
/// <param name="vA">被复制的向量</param>
/// <param name="vB">待复制的向量</param>
/// <returns>若vA和vB是等长的向量，可以复制，返回true。否则返回false</returns>
bool vec_copyto(vec vA, vec vB);

/// <summary>
/// 实现向量的复制。
/// </summary>
/// <param name="vA">待复制的向量</param>
/// <param name="vB">被复制的向量</param>
/// <returns>若vA和vB是等长的向量，可以复制，返回true。否则返回false</returns>
bool vec_copyfrom(vec vA, vec vB);

/// <summary>
/// 判断两个向量是否相等
/// </summary>
/// <param name="vA">向量A</param>
/// <param name="vB">向量B</param>
/// <param name="precision">精度要求</param>
/// <returns>如果向量每个元素的差都小于精度要求，说明两个向量相等</returns>
bool vec_isEqual(vec vA, vec vB, real_t precision);

/// <summary>
/// 判断向量是否零向量
/// </summary>
/// <param name="v">向量</param>
/// <param name="precision">精度要求</param>
/// <returns>如果向量每个元素的值都小于精度要求，说明是零向量</returns>
bool vec_isZero(vec v, real_t precision);

/// <summary>
/// 实现向量之间的加法
/// </summary>
/// <param name="A">向量A</param>
/// <param name="B">向量B</param>
/// <param name="C">向量C=A+B</param>
/// <returns>运算成功返回true。如果向量ABC的元素个数不一致，运算失败，返回false。</returns>
bool vec_add(vec A, vec B, vec C);

/// <summary>
/// 实现向量之间的减法
/// </summary>
/// <param name="A">向量A</param>
/// <param name="B">向量B</param>
/// <param name="C">向量C=A-B</param>
/// <returns>运算成功返回true。如果向量ABC的元素个数不一致，运算失败，返回false。</returns>
bool vec_sub(vec A, vec B, vec C);

/// <summary>
/// 实现向量与数的乘法
/// </summary>
/// <param name="A">向量A</param>
/// <param name="r">实数r</param>
/// <param name="B">向量B=A*r</param>
/// <returns>运算成功返回true。如果向量AB的元素个数不一致，运算失败，返回false。</returns>
bool vec_mul(vec A, real_t r, vec B);

/// <summary>
/// 计算两个向量之间的点乘。向量A和向量B的元素个数不需要一致，低维向量在高维度的分量默认为零。
/// </summary>
/// <param name="A">向量A</param>
/// <param name="B">向量B</param>
/// <returns>点乘的结果</returns>
real_t vec_dotProduct(vec A, vec B);

/// <summary>
/// 计算向量模长。
/// </summary>
/// <param name="v">向量</param>
/// <returns>模长 = sqrtf(v[0]^2+v[1]^2+v[2]^2+…+v[n]^2) </returns>
real_t vec_length(vec v);

/// <summary>
/// 计算向量模长的平方。
/// </summary>
/// <param name="v">向量</param>
/// <returns>模长平方 = (v[0]^2+v[1]^2+v[2]^2+…+v[n]^2) </returns>
real_t vec_length_square(vec v);



/// <summary>
/// 计算矢量的曼哈顿距离Manhattan distance。
/// </summary>
/// <param name="v">向量</param>
/// <returns> 曼哈顿距离 = (|v[0]|+|v[1]|+|v[2]|+…+|v[n]|) </returns>
real_t vec_manhattan(vec v);


/// <summary>
/// 计算矢量v1到v2的曼哈顿距离Manhattan distance。
/// </summary>
/// <param name="v1">向量</param>
/// <param name="v2">向量</param>
/// <returns> 曼哈顿距离 = (|v1[0]-v2[0]|+|v1[1]-v2[1]|+……) </returns>
real_t vec_manhattan_2(vec v1, vec v2);



/// <summary>
/// 随机初始化矩阵
/// </summary>
/// <param name="m">矩阵</param>
/// <param name="seed">生成随机数的种子</param>
/// <param name="min">随机数的最小值</param>
/// <param name="max">随机数的最大值</param>
void mat_randomize(mat m, int seed, real_t min, real_t max);


//实现矩阵的复制   A矩阵 => 矩阵B ，该矩阵复制给B的矩阵
bool mat_copyto(mat A, mat B);

//实现矩阵的复制   A矩阵 <= 矩阵B ，B的矩阵复制给该矩阵
bool mat_copyfrom(mat A, mat B);

/// <summary>
/// 初始化矩阵
/// </summary>
/// <param name="m">矩阵</param>
/// <param name="value">矩阵元素初始化的值</param>
void mat_init(mat m, real_t value);

/// <summary>
/// 初始化区域。左上角为坐标原点，右下方是xy的正方向
/// </summary>
/// <param name="m">矩阵</param>
/// <param name="left">区域左边（闭）</param>
/// <param name="top">区域上边（闭）</param>
/// <param name="right">区域右边（开）</param>
/// <param name="buttom">区域底边（开）</param>
/// <param name="value">初始化值</param>
void mat_initarea(mat m, int left, int top, int right, int buttom, real_t value);

/// <summary>
/// 对矩阵行矢量进行单位化
/// </summary>
/// <param name="m">矩阵</param>
void mat_row_vectorUnitization(mat m);

/// <summary>
/// 对矩阵列矢量进行单位化
/// </summary>
/// <param name="m">矩阵</param>
void mat_col_vectorUnitization(mat m);

/// <summary>
/// 单位化矩阵
/// </summary>
/// <param name="m">矩阵</param>
/// <returns>操作成功返回true。若不是方阵，操作失败，返回false</returns>
bool mat_eye(mat m);

/// <summary>
/// 实现矩阵之间的加法
/// </summary>
/// <param name="A">矩阵A</param>
/// <param name="B">矩阵B</param>
/// <param name="C">矩阵C=A+B</param>
/// <returns>运算成功返回true。如果矩阵的行列数不匹配，运算失败，返回false。</returns>
bool mat_add(mat A, mat B, mat C);

/// <summary>
/// 实现矩阵之间的减法
/// </summary>
/// <param name="A">矩阵A</param>
/// <param name="B">矩阵B</param>
/// <param name="C">矩阵C=A-B</param>
/// <returns>运算成功返回true。如果矩阵的行列数不匹配，运算失败，返回false。</returns>
bool mat_sub(mat A, mat B, mat C);

/// <summary>
/// 实现矩阵之间的乘法
/// </summary>
/// <param name="A">矩阵A</param>
/// <param name="B">矩阵B</param>
/// <param name="C">矩阵C=A*B</param>
/// <returns>运算成功返回true。如果矩阵的行列数不匹配，运算失败，返回false。</returns>
bool mat_mul(mat A, mat B, mat C);

/// <summary>
/// 实现矩阵元素之间的乘法
/// </summary>
/// <param name="A">矩阵A</param>
/// <param name="B">矩阵B</param>
/// <param name="C">矩阵C元素=A对应位置元素*B对应位置元素</param>
/// <returns>运算成功返回true。如果矩阵的行列数不匹配，运算失败，返回false。</returns>
bool mat_data_mul(mat A, mat B, mat C);

/// <summary>
/// 实现矩阵与向量的乘法
/// </summary>
/// <param name="A">矩阵A</param>
/// <param name="B">列矢B</param>
/// <param name="C">列矢C=矩阵A＊列矢B</param>
/// <returns>运算成功返回true。如果矩阵的行列数不匹配，运算失败，返回false。</returns>
bool mat_mul_vec(mat A, vec B, vec C);

/// <summary>
/// 实现矩阵转置与向量的乘法
/// </summary>
/// <param name="A">矩阵A</param>
/// <param name="B">列矢B</param>
/// <param name="C">列矢C=(A^T)*B。等效式：行矢C=行矢B*矩阵A</param>
/// <returns>运算成功返回true。如果矩阵的行列数不匹配，运算失败，返回false。</returns>
bool matT_mul_vec(mat A, vec B, vec C);

/// <summary>
/// 实现矩阵与数的乘法
/// </summary>
/// <param name="A">矩阵A</param>
/// <param name="r">实数r</param>
/// <param name="B">矩阵B</param>
/// <returns>运算成功返回true。如果矩阵的行列数不匹配，运算失败，返回false。</returns>
bool mat_mul_real(mat A, real_t r, mat B);

/// <summary>
/// 矢量乘法生成矩阵。矩阵A=列矢B＊列矢C^T
/// </summary>
/// <param name="A">矩阵A=列矢B＊列矢C^T</param>
/// <param name="B">列矢B</param>
/// <param name="C">列矢C</param>
/// <returns>运算成功返回true。如果矩阵的行列数不匹配，运算失败，返回false。</returns>
bool mat_vec_mul_vecT(mat A, vec B, vec C);



///////////////////////////////////////////////////////////////////

/// <summary>
/// 高斯消元法求逆矩阵。
/// 如果运算成功，原来的矩阵A变成单位矩阵。
/// 如果运算失败，原来的矩阵A也会被破坏。最好事先拷贝一份副本。
/// </summary>
/// <param name="A">原来的矩阵，成功求逆后变成单位矩阵</param>
/// <param name="iA">逆矩阵，运算结束后变成原来矩阵A的逆</param>
/// <param name="precision">运算精度</param>
/// <returns>如果是满秩的方阵，可以求逆，返回true。否则无法正确求逆，返回false</returns>
bool  mat_inversion(mat A, mat iA, real_t precision);

/// <summary>
/// 计算方阵的值。按照高斯消元法将矩阵分解成上三角矩阵，然后计算对角线积。
/// 如果运算成功，原来的矩阵m变成上三角矩阵。
/// 如果运算失败，原来的矩阵m会被破坏。最好事先拷贝一份副本。
/// </summary>
/// <param name="m">原来的矩阵</param>
/// <param name="precision">运算精度</param>
/// <returns>矩阵的值。如果矩阵是满秩的方阵，则可以计算值。否则返回0.</returns>
real_t mat_det(mat m, real_t precision);

/// <summary>
/// 计算矩阵的范数，即矩阵每个元素的平方和
/// </summary>
/// <param name="m">矩阵</param>
/// <returns>矩阵所有元素的平方和</returns>
real_t mat_norm(mat m);

/// <summary>
/// 实现矩阵转置
/// </summary>
/// <param name="m">矩阵</param>
/// <returns>如果矩阵是方阵，返回true，矩阵m变成转置矩阵。否则返回false，矩阵m不变</returns>
bool mat_transpose(mat m);

/// <summary>
/// 实现矩阵转置
/// </summary>
/// <param name="A">原矩阵</param>
/// <param name="AT">转置后的矩阵</param>
/// <returns>运算成功返回true。如果矩阵的行列数不匹配，运算失败，返回false。</returns>
bool mat_transpose_to(mat A, mat AT);

/// <summary>
/// 判断矩阵是否相等
/// </summary>
/// <param name="A">矩阵A</param>
/// <param name="B">矩阵B</param>
/// <param name="precision">运算精度</param>
/// <returns>A等于B，返回true。否则返回false</returns>
bool mat_isEqual(mat A, mat B, real_t precision);

/// <summary>
/// 判断该矩阵是否是上三角矩阵。只看下三角区域是否全为0，不管是不是方针。
/// </summary>
/// <param name="m">矩阵</param>
/// <param name="precision">运算精度</param>
/// <returns>m是上三角矩阵返回true，否则返回false</returns>
bool mat_isUpperTriangle(mat m, real_t precision);

/// <summary>
/// 判断该矩阵是否是下三角矩阵。只看上三角区域是否全为0，不管是不是方针。
/// </summary>
/// <param name="m">矩阵</param>
/// <param name="precision">运算精度</param>
/// <returns>m是下三角矩阵返回true，否则返回false</returns>
bool mat_isLowerTriangle(mat m, real_t precision);

/// <summary>
/// 判断该矩阵是否是对称矩阵
/// </summary>
/// <param name="m">矩阵</param>
/// <param name="precision">运算精度</param>
/// <returns>m是对称矩阵返回true，否则返回false</returns>
bool mat_isSymmetric(mat m, real_t precision);

/// <summary>
/// 判断该矩阵是否单位矩阵
/// </summary>
/// <param name="m">矩阵</param>
/// <param name="precision">运算精度</param>
/// <returns>m是单位矩阵返回true，否则返回false</returns>
bool mat_isEye(mat m, real_t precision);

/// <summary>
/// 判断该矩阵 对角线是否0或1，其他元素都为0
/// </summary>
/// <param name="m">矩阵</param>
/// <param name="precision">运算精度</param>
/// <returns>m是单位矩阵返回true，否则返回false</returns>
bool mat_isDiagonal_01(mat m, real_t precision);




///////////////////  下面算是中级内容了     ////////////////////////////


/// <summary>
/// 使用子函数的豪斯霍尔德QR分解。
/// 不要求m矩阵是方阵。A尺寸与m矩阵一样。H是和行数一样的方阵。
/// 不管A是不是满秩，分解后的H是满秩的。
/// </summary>
/// <param name="m">原矩阵=H*A</param>
/// <param name="operation_buffer">运算的缓存矩阵。m.cols行，m.rows+1列。要cols个豪斯霍尔德向量，第一个向量rows个元素。每个向量还要一个beta参数</param>
/// <param name="H">分解的正交矩阵</param>
/// <param name="A">分解的上三角矩阵</param>
/// <param name="precision">运算精度</param>
/// <returns>运算成功返回true。如果矩阵的行列数不匹配，运算失败，返回false。</returns>
bool mat_QR_Householder(mat m, mat H, mat A, real_t precision);


/// <summary>
/// 利用Givens变换实现QR分解。当该矩阵下三角区域存在大量零元的时候，Givens变换进行QR分解速度更快。
///  G是m.rows大小的方针，A的矩阵尺寸和m的矩阵一致。
/// </summary>
/// <param name="m">原矩阵=G*A</param>
/// <param name="G">分解的正交矩阵</param>
/// <param name="A">分解的上三角矩阵</param>
/// <param name="precision">运算精度</param>
/// <returns>运算成功返回true。如果矩阵的行列数不匹配，运算失败，返回false。</returns>
bool mat_QR_Givens(mat m, mat G, mat A, real_t precision);

/// <summary>
/// 使用改进“格拉姆施密特过程”的算法实现QR分解。
/// 格拉姆施密特过程不要求原矩阵是方阵，只要Q的行列数和原矩阵m一致，R是与列数相同的方阵。
/// 只有对角元全为正的时候，QR分解才是唯一的。不同方法的QR分解，分解出的Q，R不同行列向量正负号可能不同。Q的秩与原矩阵相同。
/// </summary>
/// <param name="m">原矩阵=Q*R</param>
/// <param name="Q">分解的正交矩阵</param>
/// <param name="R">分解的上三角矩阵</param>
/// <param name="precision">运算精度</param>
/// <returns>运算成功返回true。如果矩阵的行列数不匹配，运算失败，返回false。</returns>
bool mat_QR_Schmidt(mat m, mat Q, mat R, real_t precision);

/// <summary>
/// 线性方程组最小二乘法求解。m*X=mX，已知m和mX，求最佳的X。
/// 这里使用豪斯霍尔德变换进行QR分解，计算最小二乘。如果矩阵中间有零元行列，会算错。见下面代码注释。
/// </summary>
/// <param name="m">矩阵m</param>
/// <param name="operation_buffer">运算的缓存矩阵。m.rows*2+m.cols+1行，max(m.rows+1,m.cols)列。</param>
/// <param name="X">待求解的矢量X</param>
/// <param name="mX">已知的矢量mX</param>
/// <param name="is_resetX">是否重置矢量X。如果重置,计算X前会初始化为零矢量</param>
/// <param name="precision">运算精度</param>
/// <returns>操作是否成功。若矩阵不匹配，或者矩阵不是满秩，运算失败，返回false</returns>
bool mat_LinearEquation_LinearLeastSquare_Househoulder(mat m, vec X, vec mX, bool is_resetX, real_t precision);

/// <summary>
/// 线性方程组最小二乘法求解。m*X=mX，已知m和mX，求最佳的X。
/// 这里使用Schmidt变换进行QR分解，计算最小二乘。 这个算法更普适。
/// </summary>
/// <param name="m">矩阵m</param>
/// <param name="operation_buffer">运算的缓存矩阵。m.rows*2+m.cols+1行，max(m.rows+1,m.cols)列。</param>
/// <param name="X">待求解的矢量X</param>
/// <param name="mX">已知的矢量mX</param>
/// <param name="is_resetX">是否重置矢量X。如果重置,计算X前会初始化为零矢量</param>
/// <param name="precision">运算精度</param>
/// <returns>操作是否成功。若矩阵不匹配，或者矩阵不是满秩，运算失败，返回false</returns>
bool mat_LinearEquation_LinearLeastSquare(mat m, vec X, vec mX, bool is_resetX, real_t precision);



/// <summary>
/// 线性方程组最小二乘法解析求解。速度会比QR分解慢一些。
/// 如果m是横阵，这里求逆的时候，mtm的秩不够，会无法求逆。所以……横阵不要用解析法求最小二乘。
/// </summary>
/// <param name="m">矩阵m</param>
/// <param name="operation_buffer">运算的缓存矩阵。m.cols*3行，max(m.rows,m.cols)列。</param>
/// <param name="X">待求解的矢量X</param>
/// <param name="mX">已知的矢量mX</param>
/// <param name="precision">运算精度</param>
/// <returns>操作是否成功。若矩阵不匹配，或者矩阵不是满秩，运算失败，返回false</returns>
bool mat_LinearEquation_LinearLeastSquare_analytic(mat m, vec X, vec mX, real_t precision);

/// <summary>
/// 矩阵乘法最小二乘法求解。m*X=mX，已知m和mX，求最佳的X。
/// 这里将矩阵拆分成列矢进行最小二乘求解。
/// </summary>
/// <param name="m">矩阵m</param>
/// <param name="operation_buffer">运算的缓存矩阵。m.rows*2+m.cols+3 行，max(m.rows+1,m.cols)列。</param>
/// <param name="X">待求解的矩阵X</param>
/// <param name="mX">已知的矩阵mX</param>
/// <param name="is_resetX">是否重置矢量X。如果重置,计算X前会初始化为零矢量</param>
/// <param name="precision">运算精度</param>
/// <returns>操作是否成功。若矩阵不匹配，或者矩阵不是满秩，运算失败，返回false</returns>
bool mat_LinearLeastSquare(mat m, mat X, mat mX, bool is_resetX, real_t precision);


///////////  特征值分解        /////////////////


/// <summary>
/// 实现对称矩阵的本征值和本征向量的分解。 m = X*D*X^(T) 
/// </summary>
/// <param name="m">待分解的对称矩阵</param>
/// <param name="operation_buffer">运算的缓存矩阵，m.rows*3+1行，m.rows+1列。</param>
/// <param name="X">X矩阵的列向量是本征向量</param>
/// <param name="D">对角矩阵，对角元素是本征值。</param>
/// <param name="precision">运算精度</param>
/// <returns>迭代计算的次数。如果不是对称矩阵，计算失败，返回-1</returns>
int mat_EVD_Symmetric(mat m, mat X, mat D, real_t precision);

///////////  奇异值分解        /////////////////

/// <summary>
/// 实现奇异值分解
/// </summary>
/// <param name="m">m矩阵 = m_U * m_B * m_VT .</param>
/// <param name="operation_buffer">运算缓存矩阵。令：max=max(row,col)，min=min(row,col)，缓存矩阵要max*2+min*3行，max列</param>
/// <param name="m_U">m.rows*m.rows的单位正交矩阵，列矢量之间单位正交</param>
/// <param name="m_B">对角矩阵。大小和矩阵m一样，对角元素是奇异值。</param>
/// <param name="m_VT">m.cols*m.cols的单位正交矩阵，行矢量之间单位正交</param>
/// <param name="precision">运算精度</param>
/// <returns>迭代次数</returns>
int mat_SVD(mat m, mat m_U, mat m_B, mat m_VT, real_t precision);





/// <summary>
/// 测试向量相关的运算
/// </summary>
void testVector(void);


/// <summary>
/// 测试矩阵相关的运输
/// </summary>
void testMatrix(void);


void test_mat_EVD_only_Eigenvalue(void);





#endif


